Role of pinning potentials in heat transport through disordered harmonic chain 
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The role of quadratic onsite pinning potentials on determining the size (N) dependence of the 
disorder averaged steady state heat current (J), in a isotopically disordered harmonic chain con- 
nected to stochastic heat baths, is investigated. For two models of heat baths, namely white noise 
baths and Rubin's model of baths, we find that the N dependence of (J) is the same and depends 
on the number of pinning centers present in the chain. In the absence of pinning, (J)ft ~ 1/N 1 ^ 2 
while in presence of one or two pins (J) f, ~ l/N 3 ^ 2 . For a finite (n) number of pinning centers with 
2 < n <C N we provide heuristic arguments and numerical evidence to show that (J) n ~ 1/N n ~ 1 ^ 2 . 
We discuss the relevance of our results in the context of recent experiments. 
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Since the seminal paper of Anderson [l[ , the physics of 
localization in disordered systems has now been studied 
for over half a century 0, 0, 0] ■ Recently there has 
been a renewed interest in this field with a lot of work 
on some open questions such as, for example, the effect 
of interactions on localization jy, 0, HI, and the metal- 
insulator transition in two dimensions [9(. A number of 
recent experiments have also reported detailed studies 
on localization in varied systems such as heat conduction 
in a isotopically disordered nanotube (loj . electrons in 
a disordered carbon nanotube llj, photons in a waveg- 
uide and sound localization in elastic networks [l3j . 
The field is thus still filled with interesting questions and 
puzzles. Here in this paper we point out that even the 
simple problem of heat conduction in a one dimensional 
disordered harmonic lattice has surprises. 

It is well known that all the eigenstates of an elec- 
tron in a one-dimensional disordered potential are local- 
ized. The electrical current thus decays exponentially 
with wire length, making it an insulator. In contrast, 
in phononic systems, for example a disordered harmonic 
chain, long wavelength modes are extended and can con- 
duct a significant amount of heat. How good a heat con- 
ductor then is the disordered harmonic chain? The ob- 
vious question to ask is the system size (N) dependence 
of the disorder averaged steady state heat current which 
we will denote by (J), ft is expected that this has the 
form (J) ~ 1/N 1 ~ a so that the conductivity scales as 
k ~ N a . The dependence of a on choice of heat baths 
and boundary conditions has been somewhat puzzling 
and has caused some amount of confusion. We note that 
heat conduction in this system is non-diffusive and cor- 
respondingly a 7^ 0. 

We briefly review earlier work on this problem [T3 | . In 
an important work on the localization of normal modes 
in the isotopically disordered harmonic chain (IDHC), 
Matsuda and Ishii [3| (MI) showed that normal modes 
with frequencies ui ~ uid were extended. For a harmonic 
chain of length N, given the average mass m = (mi), 
the variance a 2 = ((mi — m) 2 ) and interparticle spring 



constant k, it was shown that 
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They also evaluated expressions for thermal conductivity 
of a finite disordered chain connected to (a) white noise 
baths and (b) baths modeled by semi-infinite ordered 
harmonic chains (Rubin's model of bath). In the follow- 
ing we will also consider these two models of baths and 
refer to them as model(a) and model(b). For model(a) 
MI used fixed boundary conditions (BC) and the limit of 
weak coupling to baths, while for case (b) they consid- 
ered free BC and this was treated using the Kubo for- 
malism. They found a — 1/2 in both conclusion 
which we will show is incorrect. The other two important 
theoretical papers on heat conduction in the disordered 
chain are those by Rubin and Greer [l5[ (RG)who con- 
sidered model(b) and of Casher and Lebowitz |16| (CL) 
who used model(a) for baths. RG obtained a lower bound 
(J) > l/N 1 / 2 and gave numerical evidence for an expo- 
nent a = 1/2 and this was later proved rigorously by Ver- 
heggen [lTj j. On the other hand, for model(a), CL found 
a rigorous bound (J) > 1/iV 3 / 2 and simulations by Viss- 
cher with the same baths supported the corresponding 
exponent a = — 1/2. In a more recent work [l8j], one of 
us (AD) gave a unified treatment of the problem of heat 
conduction in disordered harmonic chains connected to 
baths modeled by generalized Langevin equations and 
showed that models(a,b) were two special cases. An effi- 
cient numerical scheme was proposed and used to obtain 
the exponent a and it was established that a = —1/2 
for model(a) (with fixed BC) and a = 1/2 for model(b) 
(with free BC). It was also pointed out that in general, 
a depended on the spectral properties of the baths. 

Here we apply the same formulation as developed in 
[lH to understand in detail the role of BCs' (and more 
generally the presence of pinning potentials) on heat 
transport in the IDHC connected to either white noise 
[model(a)] or Rubin baths [model(b)]. We show that 
with the same kind of pinning, the exponent a is the 
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same for the two different bath models. The pinning po- 
tentials strongly scatter low frequency waves and hence 
can be expected to lower the heat current. Surprisingly, 
we find that even the exponent a changes with the num- 
ber of pinning centers. We also provide expressions for 
the asymptotic value of (J) for various cases. 
The Hamiltonian of the IDHC considered here is 
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where {xi,pi} denote the displacement and momentum 
of the particle at lattice site /. The random masses {mi} 
are chosen from a uniform distribution between (to — A) 
to (to + A). The strength of onsite potentials at the 
boundaries is kl . The particles at two ends are connected 
to heat baths at temperature Tl and Tr. The heat 
reservoirs are modelled by generalised Langevin equa- 
tions 18, 1|| 2(j. The steady state classical heat current 
through the chain is given by: 



J = 



k B (T L -T R ) 



dioT N (uj), (3) 



where TJv(w) 
and Z 



4tt 

4r 2 H|G 1N H| 2 , G{u) = Z- 1 /k 
[-w 2 Af + $-E(w)]/fc , 



where M and $ are respectively the mass and force ma- 
trix for the harmonic chain and G is the Green's func- 
tion of the chain connected to baths. The self-energy 
correction in the Green's function £, coming from the 
baths, is a N x N matrix whose only non-zero elements 
are En = EjVjv = £( w ) and r(w) = Jra[E]. For 
white noise baths £(w) = —ijco where 7 is the coupling 
strength with the baths, while in case of Rubin's baths 
E(w) = k{\ - muo 2 /2k - mj(to//c) 1/2 [1 - TOw 2 /(4fc)] 1/2 }. 
We have assumed that the RG bath has spring con- 
stant k and equal masses to. We note that Tn(uj) is 
the transmission coefficient of phonons through the dis- 
ordered chain. To extract the asymptotic N dependence 
of (J) we need to determine the Green's function element 
Gin(lo). It is convenient to write the matrix elements 
Z tl = -m,iu) 2 /k + l + k'/k - E/fc = - mi uj 2 /k + 2 - E' 
where £' = S/fc — k'/k + l and similarly Znn = 
mNLu 2 /k + 2 — £'. Following the techniques used in 
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we have 



\G 1N {u))\ 2 = k- 2 \A N (u)\- 2 with 



(4) 



A n (uj) = Di. N - Y/(D 2 ,n + £>i,w-i) + £' 2 £>2,at-i 

where An (us) is the determinant of Z and the matrix 
elements D^ m are given by the following product of (2 x 2) 
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where Ti 
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We note that the information about bath properties and 
boundary conditions are now contained entirely in 
while D contains the system properties. It is known that 
|-D;,m| ~ e cNu}2 for \l — m\ ~ N [3], where c is a con- 
stant, and so we need to look only at the low frequency 
(u ~ l/N 1 ' 2 ) form of E' . We now proceed to examine 
various cases. For model(a) free BC correspond to k' = 
and so E' = l — iju/k while for model(b) free boundaries 
correspond to k! — k and this gives, at low frequencies, 
E' = 1 — ^(to/A:) 1 / 2 ^. Other values of k' correspond to 
pinned boundary sites with an onsite potential k a x 2 /2 
where k Q = k' for model (a) and k Q = k' — k for model (b). 
The main difference, from the unpinned case, is that now 
i?e[E'] 7^ 1. The arguments of [l8j then immediately give 
a = 1/2 for free BC and a = -1/2 for fixed BC for both 
bath models. The arguments consisted of two parts: (i) 
it was observed numerically that the transmission coef- 
ficient at low frequencies for the ordered and disordered 
chain were almost the same, (ii) an asymptotic analysyis 
was then carried out for the ordered case, for which T 
could be obtained exactly for any bath spectral prop- 
erties (an improved version of those arguments is given 
below) . 

For the choice of parameters 7 = (to/c) 1 / 2 , the imag- 
inary part of E' is the same for both bath models, and 
we expect, for large system sizes, the actual values of the 
current to be the same in both cases. This can be seen 
in Fig. (J]) where we show the system size dependence of 
the current for the various cases. The current was evalu- 
ated numerically using Eq. ([3]) and averaging over many 
realizations (~ 4 — 100). We also show the exact asymp- 
totic forms for the current which we will discuss later. 
Note that for free BC, the exponent a = 1/2 settles to 
its asymptotic value at relatively small values (N ~ 10 3 ) 
while, with pinning, we need to examine much longer 
chains (N ~ 10 5 ). We also find that the presence of a 
single pinning centre in the IDHC is sufficient to change 
the value of a from 1/2 to -1/2 [see Fig. ©]. These 
results clearly show that, for both models(a,b), the ex- 
ponent a is the same and is controlled by the presence 
or absence of pinning in the IDHC. 

Next we try to better understand the above results. 
As mentioned before only modes ui ~ u>d are involved in 
conduction. It was noted in [l8| that in this low frequency 
regime we can approximate (T/v(w)) by the transmission 
coefficient of the ordered chain Tjy (u>) . We then obtain 
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FIG. 1: (Color online). Plot of (J) versus N for free BC 
(n = 0) and fixed BC (n — 2). Results are given for both 
models(a,b) of baths. The two straight lines correspond to 
the asymptotic expressions given in Eqs. (j 101 1 1 j) . We used 
parameters m = 1, A = 0.5, k = 1, y = 1, Ti, = 2, Tr = 1 
and k — 1. The error in the measurements is much smaller 
than the size of the symbols. 



For model(a), in the limit N — > oo is effectively given 
by [Hi: 



7<x> 2 V / Amk — m 2 uj 2 
k' 2 + (7 2 + m(k - k'))uj 2 



(7) 



We then find, for free BC (k' = 0) T°(oj) ~ 1 while for 
fixed BC (k' ^ 0), T°(w) - to 2 . Using Eq. © then 
immediately gives the asymptotic N dependence for the 
two BCs'. Our results are valid even in the weak coupling 
limit 7 < < 1 and this means that the result given by MI 
for model(a) in the weak coupling limit is incorrect. Our 
numerics supports this conclusion. We also compute the 
transmission coefficient of the ordered chain (as N — > oo) 
in the presence of a single pinning at one boundary, 



27a; 2 V Amk — m 2 u> 2 
V?4uj 4 A 2 + k>n(k>n + 4w 2 A) ' 



(8) 



with A = 7 2 + km, Q = k' — 



and we again find T° 



for k' 0. This confirms 



our numerics that the asymptotic N dependnce of (J) 
is analogous for the IDHC with single or double pinning 
centres. For model(b), the transmission coefficient of the 
ordered chain, pinned at the two boundary sites with 
k Q = k' — k is given effectively by (as N — ■> 00): 



2fc 2 sin 2 q 
2k 2 sin 2 q + k 2 ' 



(9) 



where 10 = 2{k/m) 1 / 2 sin(g/2). As expected, for k Q = 
we have T° = 1 while for k ^ 0, T° ~ uj 2 . The above 



FIG. 2: (Color online). Plot of (J) versus N for n = 1,3,4 
pinning centers for model(b). Parameters are same as in 
Fig. ([lj and for these parameters model(a) results are almost 
indistinguishable for > 10 . The straight lines have slopes 
— 1.47, —2.03 and —2.38. The error-bars shown are for disor- 
der average and are of the same order as numerical errors. 



qualitative analysis thus shows that the effect of introduc- 
ing pinning potentials is to pinch the band of conducting 
modes (between — LOd) from the zero frequency side and 
thus lower (J). 

Our asymptotic analysis also allows us to make pre- 
dictions, on the dependence of (J), on various system 
parameters such as mass variance, spring constant etc. 
Here we denote (J) Fr for (J) in the absence of pinning 
while (J) Fi represents (J) in the presence of double pin- 
nings at the boundaries. From Eq. ([6]) and the forms of 
T (u>) in various cases we get: 



(J) Fr 



A c 
A' c 



k B (T L - T R ) / 



km \ 1 / 2 



ir \Na 2 - 
, k B (T L -T R )/ km \3/2 
\Na 2 ) 



(10) 
(11) 



where c = 27(mfc) 1 / 2 /(7 2 + mk), 1 for model(a), 
model(b) respectively. For fixed boundaries we have 
d = 7(m/c) 1 / 2 jk 2 , mkfk 2 for model(a), model(b) 
respectively. A, A' are constant numbers. We 
find that for model(b) our numerical results agrees 
with an exact expression for (J) Fr due to Papan- 
icolau (apart from a factor of 27r) and this gives 
A = ir 3 / 2 f™dt [tsinh(7rf)]/[(i 2 + 1/4) 1 ' 2 cosh 2 (nt)} rj 
1.08417 (see [l7j])- We note that this differs from the 
expression given in For fixed boundaries we find nu- 
merically that A' fts 17.28 and the fit is shown in Fig. |T]). 
Based on our analytical and numerical results, we believe 
that the expression in [T7I ] is in error by a 2-7T factor. 

Till now, using numerical results and heurestic argu- 
ments, we have arrived at the result that for a IDHC, 
in the absence of any pinning potential a = 1/2 while 
the presence of one or two pinned sites changes the ex- 
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ponent to a = —1/2. This is true both for white noise 
and Rubin's bath. It is natural to now ask as to what 
happens in the presence of more number of pinning cen- 
ters. It is expected that more pinning centers will lead 
to enhanced scattering of low frequency phonons and de- 
crease the heat current but it is not obvious as to whether 
the exponent a changes. For a finite fraction of sites on 
the lattice having pinning potentials, it is known that 



(J) 



,-cN 



Here we investigate the case with a 



finite number, say n, of pinning sites. Numerically it be- 
comes difficult to determine a for n > 4 as, with more 
pins, the heat current becomes very small at large sys- 
tem sizes and numerical errors become significant. In 
Fig. ([2]) we show numerical results for n = 3,4, where 
the extra pinning potentials with k = 1 are placed in 
the bulk of the chain with equal separations. We find 
a w —1.03,-1.38 respectively for n = 3,4, which are 
clearly different from the n = 1, 2 value a = —0.5. Let 
us now see what our earlier heuristic arguments give, for 
n = 3. We again find that the low frequency behaviour of 
Ajv(w) are similar for the disordered and ordered lattices. 
Let us therefore find the form of Ajv for the ordered case. 
Let N = 2M + 1 with the 1 st , (M + l) th and N th sites 
being pinned. Except for Tm+i = T' all the other Tjs' 
are identical and given by T, say. If we denote Dm = T N 
and D' N = f M f'f M = D M f'D M , then using the fact 
that for the ordered lattice Din — sin q(N + l)/sin(g), 
where cos(q) = 1 — moj 2 /(2k), and carrying out the ma- 
trix multiplications above we find that at low frequencies 
D' 1N is larger than Din by a factor ~ l/sin(g) ~ 1/oj. 
This means that 7~n for the 3-pin case will have an extra 
factor of to 2 compared to the 2-pin case. Correspondingly 
one expects, using Eq. @, an exponent a = —3/2. The 
argument can be extended to the case of n > 2 pins (two 
of which are in the boundaries) in which case we get 



= 3/2- 



(12) 



Our numerical results for n — 3,4 (see Fig. ©) are con- 
sistent with this prediction though we are not able to 
verify the precise value of the exponent. 

Finally we note that the calculation by CL [l6| for 
the lower bound on current, in the case of two pinning 
centers (fixed boundaries) in model(a), can be extended 
to the case with more pins. The argument by CL consists 
in evaluating (D\ N ) by looking at the disorder averaged 
direct product (D ® D) = JJ^Qi) where Qi = Tj ® f t . 
In the CL case (Qi) = Q for all Q and an analysis of the 
eigenvalues of Q led to the result (D\ N ) ~ e cNu) . In our 
case, say for the case of n = 3 with an additional pinning 
at site I = M + 1, (Qu+i) — Q' is different and we have 
(D®D) = Q M Q'Q M . A careful analysis of this then 
gives (D\ N ) ~ uj- 2 e cNu}2 . Using this in Eq. © gives 
{J} 3 > Cjr °° dwuj A e- cNu>2 ~ 0[N~ 5 / 2 }, where c and C 
are constants. In general we get (J) > 0[N~ n+ ?] for 
2 < n « N. 



Quantum case: For a Hamiltonian of the form of Eq. ([2]) 
where now {xi,pi} are Heisenberg operators, the steady 
state quantum heat current through the IDHC in the 
linear response regime is given by: 



J q = 



k B {T L ~T B ) 



4tt 



dwTjv(a;)( 



2k B T> 



i2 cosech 2/ 



•2k B T' 



where T N (w) is same as Eq.© and T = (Tl + T R )/2. 
Following our derivation for the classical system we see 
that the asymptotic N dependence of (J q ) is determined 
by Tn(uj) which is here exactly the same as the clas- 
sical case. For any fixed temperature, however small, 
at suficiently large system sizes we will have HuJd << 
k B T, and hence within this cut-off frequency the factor 
(hu;/k B T) 2 cosech(huj/k B T) 2 — * 1. Hence for large sys- 
tem sizes we always get the classical result. The approach 
to the asymptotic behaviour though will be different. 

Discussion: In real experiments heat baths usually 
have a finite bandwidth making the noise correlated, as 
in Rubin's model. Here we have shown that for heat 
conduction in the IDHC these noise correlations do not 
affect the exponent a (note that a bath for which S(w) 
depends nonlinearly on ui at small frequencies can affect 
a). We have elucidated the role of boundary conditions 
and shown that the actual value of a depends on the 
number of pinned sites. Our results are also valid for 
bond disorder. We have provided explicit expressions for 
the currents which, apart from giving the system size- 
dependence, also give the dependence on various other 
parameters such as mass variance, coupling to baths etc. 
We also emphasize that heat conduction through IDHC is 
non-diffusive. Our physical understanding is as follows. 
In the presence of mass or bond disorder phonons are 
scattered coherently giving rise to localization and low 
transmission. Long wavelength phonons with to ~ tod 
[see Eq.(l)] are relatively unaffected and dominate heat 
conduction in such disordered materials. Now the in- 
troduction of pinning centers causes strong scattering of 
even the low frequency modes and, as we have shown, sig- 
nificantly reduces the current. We obtain the surprising 
and nontrivial result that the exponent a giving the sys- 
tem size dependence of current changes linearly with the 
number of pinning centers. There are now experimen- 
tal measurements of heat conduction in one-dimensional 
systems such as nanotubes and nanowires [ljj 0] and 
molecular wires 23[. At low temperatures one can ne- 



glect anharmonic effects and it will be interesting to see 
if our prediction of the strong reduction of heat current, 
by substrate potentials at localized points on a disordered 
wire, can be observed. While our results are for a sim- 
ple classical model we expect the effect of pinning to be 
quite generic and should be true for systems with more 
complicated phonon dispersions. It will be interesting to 
sec the effect of pinning potentials in heat conduction in 
two and three dimensions. 

AD thanks David Huse for useful discussions. 
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